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Аннотация 

Введение. Экологические проблемы, возникающие на мелководных водоёмах и вызываемые как природными, 
так и техногенными факторами, ежегодно наносят существенный ущерб аквасистемам и прибрежным 
территориям. Своевременно определить эти проблемы, а также пути их устранения возможно с использованием 
современных вычислительных систем. Но проведённые ранее исследования показали, что ресурсов 
вычислительных систем, использующих только центральный процессор, недостаточно для решения больших 
научных задач, в частности, по прогнозированию крупных экологических происшествий, оценке нанесенного 
ими ущерба и определению возможностей их устранения. Для этих целей предлагается использовать модели 
вычислительной системы и декомпозиции расчётной области для разработки алгоритма параллельно- 
конвейерных вычислений. Целью данной работы является создание модели параллельно-конвейерного 
вычислительного процесса для решения системы сеточных уравнений модифицированным попеременно- 
треугольным итерационным методом с использованием декомпозиции трёхмерной равномерной расчётной 
сетки, учитывающей технические характеристики используемого для расчетов оборудования. 

Материалы и методы. Разработаны математические модели вычислительной системы и расчётной сетки. 
Модель декомпозиции расчётной области выполнена с учётом характеристик гетерогенной системы. Предложен 
параллельно-конвейерный метод решения системы сеточных уравнений модифицированным попеременно- 
треугольным итерационным методом. 

Результаты исследования. На языке СОПЛА С написана программа, реализующая параллельно-конвейерный 
метод решения системы сеточных уравнений модифицированным попеременно-треугольным итерационным 
методом. Проведённые эксперименты показали, что с увеличением числа потоков время вычислений 
уменьшается и при декомпозиции расчётной сетки рациональным является разбиение на фрагменты по 


координате 2 на величину, не превышающую 10. Результаты экспериментов подтвердили эффективность 
разработанного параллельно-конвейерного метода. 

Обсуждение и заключение. По итогам проведенных исследований разработана модель параллельно- 
конвейерного вычислительного процесса на примере одного из самых трудоёмких этапов решения системы 
сеточных уравнений модифицированным попеременно-треугольным итерационным методом. Её построение 
основано на моделях декомпозиции трёхмерной равномерной расчётной сетки, учитывающей технические 
характеристики используемого в расчетах оборудования. Применение программы позволит ускорить процесс 
расчёта и равномерно по времени загрузить программные потоки. Проведенные численные эксперименты 
подтвердили математическую модель декомпозиции расчётной области. 
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Введение. В последнее время на территории Ростовской области отмечается возникновение ряда серьезных 
экологических проблем. К ним, в частности, относится эвтрофикация вод Азовского моря и Цимлянского 
водохранилища, которая вызывает рост вредоносных и токсичных видов популяций фитопланктона [1]. 
Инженерные работы в акваториях рек и морей приводят к загрязнению прилегающих территорий, изменению 
популяционной структуры биоты и ухудшению условий воспроизводства ценных и промысловых рыб. 
Изменение климата на юге России привело к увеличению количества случаев затопления некоторых территорий 
В районе Таганрогского залива и поймы реки Дон, вызванных сгонно-нагонными явлениями. В последнее 
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десятилетие в летний период несколько раз наблюдалось практически полное осушение русла реки Дон, что 
приводило к полной остановке судоходства. Чтобы спрогнозировать возникновение и развитие подобных 
случаев, спланировать пути устранения их последствий, оценить нанесенный ими ущерб, требуются 
современные программные комплексы, построенные с использованием высокоточных математических моделей, 
численных методов, алгоритмов и структур данных [2]. 

В основе математических моделей, используемых при прогнозировании природных и техногенных катастроф, 
лежат системы дифференциальных уравнений в частных производных, например уравнения Пуассона, Навье- 
Стокса, диффузии-конвекции-реакции, теплопроводности. Численное решение таких систем приводит к 
необходимости оперативного хранения больших объёмов данных (в массивах различной структуры) и решения 
систем сеточных уравнений высокой размерности, превышающих 10°. Объём оперативной памяти, требуемой 
для хранения массивов данных при численном решении только одного уравнения Пуассона для трёхмерной 
области размерностью 103х103х103 попеременно-треугольным итерационным методом, составляет более 64 Гб. 
В случае численного решения комбинированных задач требуются сотни гигабайт оперативной памяти, которые 
могут быть доступны лишь при использовании суперкомпьютерных систем. 

Проведенное ранее исследование показало, что ресурсов вычислительной системы, использующей только 
СРО, недостаточно для решения подобных научных задач [3]. Увеличение мощности и видеопамяти СРО 
позволило использовать для расчетов ресурсы видеоадаптеров [4]. Эффективность использования СРО зависит 
от применения параллельных алгоритмов для решения вычислительно-трудоемких задач водной экологии [5-7]. 
Частично решить проблемы нехватки памяти и вычислительной мощности на рабочих станциях можно 
установкой дополнительных видеоадаптеров в слоты РС1-Е Х16 непосредственно и в слоты РС-Е Х1 с помощью 
переходников РСГ-Е Х1-РС-Е Х16. Таким образом, количество видеоадаптеров, установленных на одной 
рабочей станции, можно довести до 12 [8—11]. 

Всё большую популярность в научном сообществе приобретают гетерогенные вычислительные системы, 
которые позволяют использовать ресурсы СРИ и СРИ совместно. Применение таких систем дает возможность 
уменьшить время расчета научных задач [12—14]. Однако эффективное использование гетерогенной 
вычислительной среды предполагает модернизацию математических моделей, алгоритмов и программ, их 
численно реализующих. Гетерогенная система позволяет организовать процесс вычислений в параллельном 
режиме. При этом должны быть учтены принципиальные различия в построении программных систем, 
использующих совместно СРО и СРО. 

Материалы и методы. Опишем предложенные математические модели вычислительной системы, расчетной 
сетки, а также метод декомпозиции расчетной области. 


Пусть р — множество технических характеристик вычислительной системы, тогда: 
Р=рц)р?() Вз, (1) 
где Ш! — подмножество характеристик центральных процессоров (СРО) вычислительной системы; 


р? — подмножество характеристик видеоадаптеров (СРО) вычислительной системы; ОЗ — подмножество 
характеристик оперативной памяти. 


Ш = (=, 412, 4:3, 41“), (2) 
где 4" — суммарное число ядер СРО; 41.2 — количество потоков, одновременно обрабатываемых одним ядром 


процессора СРО; 413 — тактовая частота, МГц; 4“ — частота шины центрального процессора, МГц. 


22= 0 2? = (4? | ЗАсро © Кори» 4? =1,„ } ? ы 


Асьи 


АСРИ ЕКбьу 


где Ксы, ={1,..,М ри} — множество индексов видеоадаптеров; №», — количество видеоадаптеров 
вычислительной системы; Ас», — индекс видеоадаптера. 


Каждый видеоадаптер представим в виде кортежа: 


В = Е > Ч ), (4) 
где о — объём видеопамяти видеоадаптера с индексом Ас», Гб; - — количество потоковых 
мультипроцессоров. 

3 = (аз, 432, (5) 
где 43' — суммарный объём оперативной памяти, Гб; 43? — тактовая частота оперативной памяти, МГц. 
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Пусть 5 — множество программных потоков, задействованных в вычислительном процессе, тогда: 


5=51 152, 

К и (6) 

я Ц 5. 5, о. 
где 5! — подмножество программных потоков, реализующих процесс расчета на СРО; 52 — подмножество 
потоковых блоков СОРА, реализующих процесс расчета на потоковых мультипроцессорах СРИ; 
№ — количество задействованных программных потоков СРИ; 5, — подмножество потоковых блоков 


СОРА, реализующих процесс расчета на потоковых мультипроцессорах СРО с индексом А 
К 


№,  —_ количество задействованных потоковых блоков СОРА, реализующих процесс расчета на потоковых 


АСРИ 


мультипроцессорах СРО с индексом Аср. 


СРИ › 


о Мы — множество индексов СРО; М№;„, — количество СРО в вычислительной системе; 


Пусть Е — множество идентификаторов программных потоков. Тогда для идентификации программных 
потоков в вычислительной системе каждому элементу множества программных потоков 3 поставим в 


соответствие кортеж е из двух элементов: 
У5Еб ееЕ: е=(и,,п,), (7) 


где п, — индекс вычислительного устройства в вычислительной системе; п, — индекс программного потока 


СРО или потокового блока СРО. 


ны | 5Е5! (8) 
а—_ 2 
Ксьуь 5 Е 5? ? 
_]|Кя› 5е8 9 
п, 1 К „о з\Е 52 . ( ) 
Зри РО 
Возьмём расчётную область со следующими параметрами: /, — характерный размер по оси Ох; [, — по 
оси Оу; | — пооси 02. Сопоставим с указанной областью равномерную расчётную сетку следующего вида: 
И =х =, у, = №,,2, =, 
1=0,п, 1 =0,и, 1, =0,п, —1; (10) 
(и, 1), = 1, (и, 1, =1,,(, 1), =1,}, 
где Л, й,, Й. — шаги расчётной сетки по соответствующим пространственным направлениям; и,, и, и, — 
количество узлов расчётной сетки по соответствующим пространственным направлениям. 
Тогда множество узлов расчётной сетки представим в виде: 
С = вый =0,п, —1,/ =0,7, —1^ =0, п. 1, (11) 
“лк = (1%. У 2%), 
где <, ,‚, — Узел расчётной сетки. 
Число узлов расчётной сетки №; вычисляется по формуле: 
М =п,-п,-п,. (12) 
Под подразделом расчётной сетки С“ с С (далее — подраздел) будем понимать подмножество узлов 
расчётной сетки С. 
= Ц 6% ={в* [Эк еК,, 84 еС,], Г] 04 =©, (13) 


®ЕКи ЕКи 
где К, = {1,..,М№ _ — множество индексов подразделов С* расчётной сетки С; № — количество подразделов С\; 
К к» № в > №; № — множество натуральных чисел; К, — индекс подраздела С\. 
Так как С“ < С ‚тогда: 
@® = а =0,м, 1.1 =0, ий БА =б,м, 1 (14) 


а 


где Яр —_ узел подраздела К, ‚ знак => обозначает принадлежность к подразделу; 7 —_ индекс узла подраздела 
м 


К, по координате у; п, 


—Й_ количество узлов в подразделе К, по координате у. 


а 


и. (15) 
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где п) 


— количество узлов по координате у В -го раздела. 
Под блоком расчётной сетки С*“ю (далее — блок) будем понимать подмножество узлов расчётной сетки 


подраздела С“. 


Ге = |9) С.ь = {ди ЭА, = К, > ды = 6%, }, П Сн.ь = [62 р (16) 
ЕК, 1 го ЕК 
1,2 1,42 
где К, = ть ,, Е — множество индексов блоков С“ подраздела (*; №, ‚, — количество блоков С**® ; 


Кик, М, СМ; А, — индекс блока С“ подраздела С*. 


Так как С*.% < С\, тогда: 


= В.Ю :_ Е и, _ ви — 
еси Е 2 =0,п, 1,7 =0, пу” —1,К =0,и, 1, (17) 
где 2“ — узел блока К,К,; знак ^ обозначает принадлежность к блоку; 7 — индекс узла блока АК, по 


К 


РА 


м, 


координате у; п, —Щ_ количество узлов в блоке К, по координате у. 


й =, У» С 
А ^ 18 
х =, у= хх п\ + } -Й,ь 2, ЁЙ,, (8) 
Вы = ^ 


В ,Б> 
у 


Под фрагментом расчётной сетки С“. (далее — фрагмент) будем понимать подмножество узлов расчётной 


где п›”” — количество узлов блока В.,Б,. 


сетки блока С“. подраздела С“. 


Са = | у Сь.вь = [анны |3, е К, ‚ав © быв}, 
П 6% = в, (19) 
БЕК 
тде А ьь= тым. т — множество индексов фрагментов С**^% блока С“ подраздела С*; 
№ к — Количество фрагментов С\ >; К, Мк СМ; А — индекс фрагмента С“** блока С\® 


подраздела С\. 
Каждому индексу А, фрагмента С*“*%* поставим в соответствие кортеж индексов 1 К; ь предназначенный 


для хранения координат фрагмента в плоскости хО2 ‚› где К. — индекс фрагмента по координате х; К; — индекс 


фрагмента по координате 2. 


В (20) 
где К, — индекс фрагмента по координате х; А, — индекс фрагмента по координате 2. 
Количество фрагментов Сл.» блока С. вычислим по формуле: 
К, =К, 'Ке ь (21) 
где К, — количество фрагментов по оси Ох; К, — количество фрагментов по координате 2. 


Так как С“ < С“® ‚ тогда: 


Съыь = в т=0,й, -1 71=0,й,- #=0,И, —1} : (2) 
где 5, ‚Е — Узел фрагмента; знак _. обозначает принадлежность к фрагменту; 1, К — индексы узла фрагмента 
по координатам х, 2; й,, И, — количество узлов расчётной сетки в фрагменте по координатам х, 2; [,, 
1. — размеры фрагмента по координатам х, 2. 

Е: и (х; У,» а: 
ИО А ИО (23) 
х, = руль + Й» У, = ЛЩ,, 2, = 2. ль +Ё й., 


где и, — количество узлов В -го фрагмента. 


Введем множество сопоставлений блоков расчётной сетки программным потокам М! 


М'= 0 им! , (24) 


К, 
К ЕКА \ ЕК, 


где М; „ — элемент множества М'. 
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1 К а 
Пусть М, „, — сопоставление блока С^ программному потоку 5, „, ‚ тогда: 
1 — А, № 
М! = (644,8, ), (25) 
где 5, ‚ е5 — программный поток, вычисляющий блок С\®. 


В процессе решения задач гидродинамики на трехмерных расчетных сетках большой размерности 
необходимы высокопроизводительные вычислительные системы и огромные объемы памяти для хранения 
данных. Ресурсов одного вычислительного устройства недостаточно для вычислений и хранения трехмерной 
расчетной сетки со всеми ее данными. Для решения этой проблемы предложены различные способы 
декомпозиции расчетных сеток с последующим применением параллельных алгоритмов расчета в гетерогенных 
вычислительных средах [15]. 

Для декомпозиции расчётной сетки необходимо учитывать производительность вычислительных устройств, 
участвующих в расчётах. Под производительностью будем понимать количество узлов расчётной сетки, 
рассчитываемых с помощью заданного алгоритма в единицу времени. 

Предположим, что все вычислительные устройства используются для расчётов. Тогда суммарная 
производительность вычислительной системы Р. вычисляется по формуле: 


Мсви 
Б; =Ери Му + > РМ, (26) 
Ь 


где Рори —_ производительность одного потока СРО; М —_ число программных потоков, реализующих процесс 


расчета на СРО; и — производительность СРО с индексом Б на одном потоковом мультипроцессоре; 


№, — количество потоковых блоков СОПА, реализующих процесс расчета на потоковых мультипроцессорах СРО. 


Тогда рассчитать количество узлов расчётной сетки п в подразделе по координате у для каждого СОРИ с 


индексом Б можно по формуле: 


Р? 
= | би |, 
те |. 7) 
х 
В процессе вычисления по формуле (27) получим остаток — некоторое количество узлов расчётной сетки. 


Эти узлы будут располагаться в оперативной памяти. Количество оставшихся узлов ив по координате у 


рассчитывается по формуле: 


М ри 
п =т- У. (28) 


= 
Для вычисления количества узлов по координате у в блоках расчётной сетки, обрабатываемых ПОТОКОВЫМИ 


мультипроцессорами СРО, воспользуемся формулами: 


Пт | т 
ПЕ. = (29) 


где и’, — количество узлов по координате у в блоках расчётной сетки, обрабатываемых потоковыми 


мультипроцессорами СОРИ с индексом Ь › кроме последнего блока; пт —_ количество узлов по координате у 


в последнем блоке расчётной сетки, обрабатываемом потоковыми мультипроцессорами СОРИ с индексом Ь. 
Для вычисления количества узлов расчётной сетки По координате у в блоках, обрабатываемых 
программными потоками, реализующими процесс расчета на СРО, воспользуемся формулами: 
ИсРИ 
Ва , 


Пт = 
УСТ М —1 (30) 


Е 
— ИСРИ — : _ 
Пуст = ПУ Пуст (м ы 1), 


где п: — количество узлов расчётной сетки по координате у, обрабатываемых программными потоками СРО, 


УСТ 


кроме последнего потока; ип —_ количество узлов расчётной сетки по координате у, обрабатываемых 


УСТЕ 
программными потоками СРО, в последнем потоке. 
Рассчитаем количество фрагментов расчётной сетки по координате у: 


М ри 
мЕМ + №. (31) 


9 


Пусть задано количество фрагментов №! и М№/ по координатам х и 2 соответственно. Тогда количество 


узлов расчётной сетки по координате х вычисляется по формулам: 
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№/-—1 | (32) 
прп, и! (М1) 


пу = 


где п/ —Ш_ количество узлов расчётной сетки по координате х во всех фрагментах, кроме последнего; 


пл. — количество узлов расчётной сетки по координате х в последнем фрагменте. 


Аналогично вычисляется количество узлов расчётной сетки по координате 2: 
М1 (33) 


где и/ — количество узлов расчётной сетки по координате2 во всех фрагментах, кроме последнего; 
пл. — количество узлов расчётной сетки по координате 2 в последнем фрагменте. 
Опишем модель параллельно-конвейерного метода. Пусть на М! необходимо организовать параллельный 


процесс вычислений некоторой функции ЕЁ, причем вычисления в каждом фрагменте С“? зависят от значений 


в соседних фрагментах, каждый из которых имеет хотя бы один из индексов по координатам х, у и 2 на 


единицу меньший, чем у текущего (рис. 1). 
Для организации параллельно-конвейерного метода введём множество кортежей А, задающих 


соответствия а между идентификаторами программных потоков е, обрабатывающих фрагменты С’*.ю.в , 


номерам шагов параллельно-конвейерного метода и: 


УеЕЕЁ ЗаеЛ: а= (е, би, т), (34) 
где г=1, №, — номер шага параллельно-конвейерного метода; №, — число шагов параллельно-конвейерного 
метода, вычисляемое по формуле: 

М, = ММ + МГ -1. (35) 


Полная загрузка всех вычислителей в предлагаемом параллельно-конвейерном методе начинается с шага 


Поозтлкг = №/ и заканчивается на шаге Идзтор = М/М№/ . При этом общее количество шагов с полной загрузкой 


вычислителей М№„.„ составит: 
Мрак = Подтор — Поовтакг +1= МУ МУ — Му +1. (36) 
Время вычислений некоторой функции Ёпараллельно-конвейерным методом запишем в виде: 
т. = х тах(Т.) , (37) 
где Т, — вектор значений затрат времени на обработку фрагментов в параллельном режиме. 
0 1 2 3 №1 
ео — 0 '=0 7=1 7=2 '—3 й 


(№, 1 › 1 ›0) 


3,1,0 
=М 


(0,1,0) | (1,1,0) | (21,0) 
И 2 | 3 

(0,2,0) | (1,2,0) | (2,2,0) | (3,2,0) (М№,-1,2,0) 
в > — | р 

(0,3,0) | (1,3,0) | (2,3,0) (№,-1,3,0) 
3 ——> 3 И ›= М2 


Рис. 1. Параллельно-конвейерный вычислительный процесс 
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Результаты исследования. Вычислительные эксперименты проводились на высокопроизводительной 
вычислительной системе К-60 Института прикладной математики им. М. В. Келдыша Российской академии 
наук. Использовалась секция с СРО, каждый узел которой оснащён двумя процессорами ше Хеоп Со!4 6142 у4, 
четырьмя видеоадаптерами №1 а Уоа СУ100ОТ, и 768 Гб оперативной памяти. 

Вычислительный эксперимент состоял из двух этапов —Щ подготовительного и основного. На 
подготовительном этапе проверялась корректность декомпозиции расчетной области на подразделы, блоки и 
фрагменты путем поэлементного сравнения значений в узлах исходной сетки и в фрагментах, полученных в 
результате декомпозиции. Затем проверялась работа алгоритма управления потоками, в процессе которого 
фиксировалось время, затрачиваемое на расчет 1, 8, 16 и 32 фрагментов расчетной сетки размерностью 50 узлов 
по пространственным координатам х, у и 2 тем же количеством потоков СРИ №, итерационным 


попеременно-треугольным методом в параллельном режиме. Выполнялось 10 повторов с вычислением среднего 
арифметического Т, и стандартного отклонения с. По полученным данным вычислялось время Т'! =Т, / Ма, 
затрачиваемое каждым потоком на обработку одного фрагмента расчетной сетки и ускорение ЕЁ = Т\(М№,)/Ти(О, 
равное отношению времени 7'(М№„) обработки одного фрагмента №. потоками к соответствующему времени 
обработки одним потоком 71(1). Экспериментальные данные приведены в таблице 1. Проведённый эксперимент 


показал, что стандартное отклонение имеет наименьшее значение в случае использования 32 параллельных 
потоков СРИ и составляет 0,026 мс, то есть использование 32 параллельных потоков СРИ при расчёте 
32 фрагментов расчётной сетки даёт более равномерную по времени загруженность программных потоков, что в 
целом повышает эффективность работы вычислительного узла. При этом среднее значение расчета одного 
4,14 мс. 
линейная Е =0,603+0,804/№., с коэффициентом детерминации, равным 0,99. Получили, что с увеличением 


фр агмента составило Зависимость ускорения Е от числа потоков получилась 


числа потоков ускорение разработанного алгоритма возрастает. Это свидетельствует об эффективном 
использовании подсистемы при работе с памятью. 


Таблица 1 
Результаты подготовительного этапа вычислительного эксперимента 
Ма тлах (Т,) , мс о, мс Т,, ме Е 
1 3,38 0,141 3,38 1,00 
8 3,66 0,042 0,46 7,39 
16 3,94 0,028 0,25 13,73 
32 4,14 0,026 0,13 26,13 


На основном этапе вычислительного эксперимента трёхмерная расчетная область, имеющая размеры 1600, 
1600, 200 по пространственным координатам х, у и 2 соответственно, разбивалась на 32 фрагмента по 50 узлов 


по каждой из координат хи у. Разбиение на фрагменты по координате 2 приведено в таблице 2. Для каждого 


варианта декомпозиции с десятикратной повторностью замерялось время обработки всей расчетной сетки 
предлагаемым параллельно-конвейерным методом и вычислялось его среднее значение я Ускорение В 


вычислялось как отношение Г’ ко времени Т расчета последовательной версией алгоритма, равному 6963 мс. 


рт эт 


Получено уравнение регрессии Е», = 7,35+1,97 1(№/) с коэффициентом детерминации, равным 0,94. Анализ 


результатов основного этапа вычислительного эксперимента показал существенное замедление роста Е при 


рт 
М№/ > 10. Поэтому делаем вывод о том, что оптимальным является разбиение на фрагменты по координате 2 на 
величину, не превышающую 10. 
Таблица 2 
Результаты основного этапа вычислительного эксперимента 


№/ п/ Ти, мс р 
1 200 1033,20 6,74 
2 100 779,00 8,94 
4 50 651,90 10,68 
8 25 588,35 11,84 

20 10 550,22 12,66 
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Обсуждение и заключение. В результате проведённых исследований разработана модель параллельно- 
конвейерного вычислительного процесса на примере одного из самых трудоёмких этапов решения системы 
сеточных уравнений модифицированным попеременно-треугольным итерационным методом. Её построение 
основано на моделях декомпозиции трёхмерной равномерной расчётной сетки, учитывающей технические 
характеристики используемого в расчетах оборудования. 

Результаты, полученные в ходе вычислительных экспериментов, подтверждают эффективность 
разработанного метода. Подтверждена и корректность декомпозиции расчетной области на подразделы, блоки и 
фрагменты. Проверена работа алгоритма управления потоками, при этом выявлено, что стандартное отклонение 
имеет наименьшее значение в случае использования 32 параллельных потоков СРО и составляет 0,026 мс, то есть 
использование 32 параллельных потоков СРУ при расчёте 32 фрагментов расчётной сетки даёт более 
равномерную по времени загруженность программных потоков. При этом среднее значение расчета одного 
фрагмента составило 4,14 мс. 

Результаты обработки замеров времени расчетов предлагаемым параллельно-конвейерным методом показали 
существенное замедление роста ускорения при разбиении на фрагменты по координате 2 при М№/ > 10. Получили, 


что оптимальным является разбиение на фрагменты по координате 2 на величину, не превышающую 10. 
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